function f = GetOccupation(E, T)
    f = (exp(E/(1.380649e-23*T)) + 1).^(-1);
end